function hfile = srf_AddMarkers(hfile, crd, opts)
% SRF::AddMarkers  - add arbitrary markers at coordinates
%
% FORMAT:       [srf] = srf.AddMarkers(crd [, opts]);
%
% Input fields:
%
%       crd         Cx3 coordinates
%       opts        struct with optional fields
%        .color     1x3 or Cx3 RGB color, otherwise red
%        .shape     one of {'sphere'}
%        .size      1x1 or Cx1 doulbe, default 3 units
%        .tal2bvc   flag, apply TAL to BV transform (default: true)
%
% Output fields:
%
%       srf         altered object

% Version:  v0.7f
% Build:    9021300
% Date:     Feb-13 2009, 12:02 AM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% load unit sphere
global bvqxfile_config;
if ~isfield(bvqxfile_config.loadedobjs, 'srf_AddSurfacePoints')
    try
        ufn = [bvqxtools_path('srf') '/unitsphere320.srf'];
        [usrf, usrfl] = BVQXfile(ufn);
        usrfbc = bvqxfile_getcont(usrf.L);
        if ~isBVQXfile(usrf, 'srf')
            error('BAD_SURFACE');
        end
        bvqxfile_config.loadedobjs.srf_AddSurfacePoints = struct;
        bvqxfile_config.loadedobjs.srf_AddSurfacePoints.unitsphere = usrfbc;
        if usrfl
            bvqxfile_clear(usrf.L);
        end
    catch
        error( ...
            'BVQXfile:FileNotFound', ...
            'Required unit sphere surface not found: ''%s''.', ...
            ufn ...
        );
    end
else
    try
        usrfbc = bvqxfile_config.loadedobjs.srf_AddSurfacePoints.unitsphere;
    catch
        error( ...
            'BVQXfile:InternalError', ...
            'Unit sphere surface removed from global config.' ...
        );
    end
end

% check arguments
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf') || ...
   ~isa(crd, 'double') || ...
    size(crd, 2) ~= 3 || ...
    any(isinf(crd(:)) | isnan(crd(:)) | crd(:) < -128 | crd(:) > 256)
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
srfbc = bvqxfile_getcont(hfile.L);
ncrd = size(crd, 1);
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'color') || ...
   ~isa(opts.color, 'double') || ...
    size(opts.color, 2) ~= 3 || ...
   ~any(size(opts.color, 1) == [1, ncrd]) || ...
    any(isinf(opts.color(:)) | isnan(opts.color(:)) | opts.color(:) < 0)
    opts.color = [255, 0, 0];
else
    opts.color = min(255, opts.color(:, 1:3, 1));
    if all(opts.color(:) <= 1)
        opts.color = opts.color * 255;
    end
    opts.color = round(opts.color);
end
if size(opts.color, 1) == 1
    opts.color = repmat(opts.color, [ncrd, 1]);
end
if ~isfield(opts, 'size') || ...
   ~isa(opts.size, 'double') || ...
    numel(opts.size) ~= length(opts.size) || ...
   ~any(numel(opts.size) == [1, ncrd]) || ...
    any(isinf(opts.size) | isnan(opts.size) | opts.size <= 0)
    opts.size = 3;
else
    opts.size = opts.size(:);
    opts.size(opts.size > 16) = 16;
end
if numel(opts.size) == 1
    opts.size = repmat(opts.size, [ncrd, 1]);
end
if ~isfield(opts, 'tal2bvc') || ...
   ~islogical(opts.tal2bvc) || ...
    numel(opts.tal2bvc) ~= 1
    opts.tal2bvc = true;
end
if opts.tal2bvc
    crd = bvcoordconv(crd, 'tal2bvc', aft_BoundingBox(hfile));
end

% set to RGB color mode and make copy of unit sphere coordinates
usrfbc.VertexColor(:, 1) = NaN;
usrfc = usrfbc.VertexCoordinate;
nvert = size(usrfc, 1);
usrfc = usrfc - repmat(usrfbc.MeshCenter, [nvert, 1]);

% make new srf
nsrf = BVQXfile('new:srf');
nsrfbc = bvqxfile_getcont(nsrf.L);
nsrfbc.NrOfVertices = 0;
nsrfbc.NrOfTriangles = 0;
nsrfbc.MeshCenter = srfbc.MeshCenter;
nsrfbc.VertexCoordinate = zeros(0, 3);
nsrfbc.VertexNormal = zeros(0, 3);
nsrfbc.VertexColor = zeros(0, 4);
nsrfbc.Neighbors = cell(0, 2);
nsrfbc.TriangleVertex = zeros(0, 3);
nsrfbc.NrOfTriangleStrips = 0;
nsrfbc.TriangleStripSequence = zeros(0, 1);

% get coordinate table of fiducials
fidt = [crd, opts.size, opts.color];

% first add all spheres to new srf
for sc = 1:ncrd
    fid = fidt(sc, :);
    fidc = repmat(fid(1:3), [nvert, 1]);
    usrfbc.VertexCoordinate = fid(4) * usrfc + fidc;
    usrfbc.VertexColor(:, 2:4) = repmat(round(fid(5:7)), [nvert, 1]);
    cnv = size(nsrfbc.VertexCoordinate, 1);
    nsrfbc.VertexCoordinate = ...
        [nsrfbc.VertexCoordinate; usrfbc.VertexCoordinate];
    nsrfbc.VertexNormal = ...
        [nsrfbc.VertexNormal; usrfbc.VertexNormal];
    nsrfbc.VertexColor = ...
        [nsrfbc.VertexColor; usrfbc.VertexColor];
    nnei = ...
        [nsrfbc.Neighbors; usrfbc.Neighbors];
    nsrfbc.TriangleVertex = ...
        [nsrfbc.TriangleVertex; (usrfbc.TriangleVertex + cnv)];
    for vc = (cnv+1):size(nnei, 1)
        nnei{vc, 2} = nnei{vc, 2} + cnv;
    end
    nsrfbc.Neighbors = nnei;
end

% make nsrf settings
nsrfbc.NrOfVertices = size(nsrfbc.VertexNormal, 1);
nsrfbc.NrOfTriangles = size(nsrfbc.TriangleVertex, 1);
bvqxfile_setcont(nsrf.L, nsrfbc);

% add spheres to srf
nfile = srf_Combine(hfile, nsrf, struct('type', 'wholebrain'));
bvqxfile_setcont(hfile.L, bvqxfile_getcont(nfile.L));
bvqxfile_clear(nfile.L);
